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Abstract 


O?  fird  H  {S  ' • 


N  points  in  a  square  of  area  A  may  be  sorted  according  their  images  under 
a  spacefilling  mapping  to  give  a  tour  of  length  at  most  If  the  points 

are  statistically  independent  under  a  smooth  distribution,  with  N  large,  then 
the  tour  will  be  roughly  25%  longer  than  optimum  (and  a  simple  enhancement 
reduces  this  to  15%)*  The  algorithm  is  easily  coded:  a  complete  BASIC  program 
is  included  in  the  appendix.  Since  the  algorithm  consists  essentially  of 
sorting,  points  are  easily  added  or  removed.  Our  method  may  also  be  used 

...  _ — ..  -  -  - — * 
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1.  Introduction 

The  travelling  salesman  problem  (TSP)  is  to  construct  a  circuit  of  minimum 
total  length  that  visits  each  of  N  given  points.  Even  In  the  plane,  this  problem 
NP-complete  (5].  Karp  [4J  has  given  an  asymptotically  efficient  heuristic,  but 
it  Is  difficult  to  code  and  its  effort  has  a  large  constant  factor.  Bentley  and 
Saxe  [2]  gave  an  efficient  implementation  of  the  nearest  neighbor  heuristic, 
but  it  requires  a  special  data  structure.  We  give  a  faster,  simpler  heuristic 
that  performs  comparably. 

2.  Overview  of  the  Method 

Let 

C  -  {0  (  0  <  0  <  1}  (1) 

denote  the  unit  circle,  so  that  0  e  C  represents  a  point  on  the  circle  0 
revolutions  removed,  clockwise,  from  a  fixed  reference  point.  Also  let 

S  ■  {(x,y)  |0<*<1,  0<y<ll  (2) 

denote  the  unit  square,  and  suppose  that  we  are  given  a  continuous  mapping  P 
from  C  onto  S.  Such  mappings  were  first  constructed  by  Peano  and  Hilbert  in 
the  1890's  and  are  known  as  "spacefilling  curves."  (See,  e.g.,  Hobson  [3, 
pp.  451-458], 

Suppose,  moreover,  that  lim^^pOS)  *  p(0).  Then,  as  0  ranges  from  0  to  1, 
p(6)  traces  out  a  "tour"  of  all  the  points  in  S.  Given  N  points  in  S  to  be 
visited,  a  reasonable  strategy  is  to  sequence  them  as  they  appear  along  the  space' 
filling  curve.  In  short,  we  sequence  them  according  to  their  inverse  image 
under  p. 
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The  specific  fora  of  ^  is  not  essential  to  our  presentation  at  this 
time,  and  is  deferred  to  Section  3.  He  need  only  two  properties: 

(PI)  An  inverse  iaage  of  can  be  easily  computed.  Specifically, 
if  x  and  y  have  k-bit  binary  representations,  then  a  9  satisfying  (x,y)  “  <»(6) 
aay  be  computed  in  0(k)  operations.  Note  that  although  many  such  9  nay  exist, 
we  must  compute  only  one. 

(P2)  There  is  a  concave  function  f(*)  on  (0,1),  with  f(0)  *  0  and 
f (A)  *  f(l-A),  such  that 

I |  *(e>  -  <K0')  | |  <.  f(|e-e'|)  (3) 

where  |  j  •  1 1  denotes  the  metric  on  S  with  respect  to  which  the  tour  is  to 
be  minimized. 

3.  Routing  Problems  in  the  Circle 

The  basic  idea  in  our  approach  is  to  solve  routing  problems  in  C  rather 
than  S,  taking  as  the  distance  between  two  points  in  C  the  upper  bound  f([9-9'|) 
on  the  distance  between  their  images  under  d*.  The  problem  in  C  has  a  great 
deal  of  structure,  attributable  to  the  concavity  of  f: 

Proposition  1:  (Triangle  Inequality)  Let  0^  —  ®2  ~  03*  ^en  +  ^03-02^ 

>  f (0^-01) . 

Proposition  2:  (Crossing  Elimination)  Let  0^  <_  ©2  <.  9^  £  0^.  Then 
(i)  f (0^-Oj)  +  f(04-02)  >  f(02-@1)  +  f(04>e3) 

Ui)  f(e3-e1)  +  f(04-02>  >  f(e3“e2)  +  f(Vei> 

1 

Proofs:  We  appeal  to  the  following  standard  inequality  for  concave  functions: 

f(a)  +  f(d)  £  f (b)  +  f(c)  if  a  £  b,c;  b,c  <  d;  a  +  d  ■  b  +  c. 

To  prove  Proposition  1,  let  a  ■  0,  b  ■  92  -  0^,  c  ■  63  -  02»  and  d  •  0_  -  0^. 


To  prove  Proposition  2  (i),  let  a  “  ®2  "  ®1*  P  “  ~  ei’  c  *  1  ~  ”  ®2^’ 

d  **  1  -  (8^  -  8^),  and  recall  that  f(A)  -  f(l-A).  The  proof  of  part  (ii)  is 
similar.  □ 

Remark:  Our  metric  is  more  general  than  | 6-8 ' | ,  and  the  triangle  inequality 
will  usually  be  a  strict  inequality! 

It  follows  from  these  propositions  that  an  optimal  tour  on  C,  under 
metric  f ,  is  obtained  by  visiting  the  points  in  a  sequence  from  smallest  6 
to  largest  8. 

Variations  of  the  TSP  may  also  be  easily  solved  on  C.  For  example.  If 

.  *  ***** 
all  pointa  are  to  be  visited  starting  at  8  and  ending  at  8  ,6  <8  ,  then 

we  know  that  all  points  between  8*  and  8  **  must  be  visited  in  an  increasing 

sequence,  the  remaining  points  must  be  visited  in  decreasing  sequence  from 
*  ** 

6  to  0  and  then  from  1  to  8  .  The  optimal  lnterleaflng  of  these  two  sequences 

2 

is  obtained  in  Q(N  )  operations,  by  dynamic  programing.  The  TSP  path  problem 
with  one  free  endpoint  is  similarly  solvable. 

Note  that,  once  the  optimal  tour  has  been  obtained  for  any  set  of  points 
on  C,  it  has  been  obtained  for  any  subset  (not  necessarily  consecutive)  of  those 
points*  ' 


4.  The  Spacefilling  Curve 

The  approach  outlined  in  Section  2  la  valid  for  any  spacefilling  curve 
f  satisfying  (PI)  and  (P2) .  He  now  describe  one  such  curve  for  which  the 

I 

bound  (3)  is  particularly  tight.  This  curve  is  recursively  defined  by  dividing 
S  into  four  identical  suhsquares,  constructing  a  curve  that  fills  each  sub¬ 
square,  and  joining  them  at  the  center  of  S  (see  Figure  1).  The  Intuition 
behind  this  definition  is  to  join  a  point  with  its  imsdiate  neighbors  before 
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proceeding  to  a  new  region.  In  this  respect  It  resembles  the  partitioning 
algorithm  of  Karp  [4]. 

•  In  order  to  translate  this  idea  into  a  mathematical  equation,  we  first 
number  the  vertices  of  S: 


q0  *  (0,0),  qx  =  (0,1),  q2  -  (1,1),  q.,  -  (1.0) 
and  adopt  the  convention: 


(4) 


ip(i/ 4)  <*  i  ■  0,1, 2, 3 


(5) 


Accordingly,  4»(G)  will  cover  the  subsquare  of  S  whose  outside  vertex  is 
when  | e— i / 4 |  1/8  (or  >_  7/8  when  i  -  0) .  Each  subsquare  covering  must  be 

rotated  so  chat  it.  begins  and  ends  at  the  center  of  S.  Consequently  we 
may  express  ip  recursively  as  the  solution  of 


V*  CO) 


*sli|»(fract(4e  +  (6-i)/4))  + 


i 


int(49  +  *s)  mod  4 


(6) 


where  the  term  scales  each  subsquares,  thfe  argument  of  f ract ( • )  reorients  it, 
and  +q^  translates  it  to  the  desired  position  in  S. 

If  we  view  (6)  as  a  fixed  point  identity,  ip  Tq» ,  then  T  is  seen  to  be  a 

contraction  operator,  and  so  there  is  a  unique  function  y  satisfying  (6).  More¬ 
over  r’V  •+  p  for  any  Initial  approximation  if  of  p.  A  sequence  of  approximations 
starting  from  (5)  is  shown  in  Figure  2. 

The  same  argument  may  be  carried  out  in  reverse  to  show  that  there  is  a 
function  ft  8  *  C  such  that  i|>($(x,y))  *  (x,y).  This  function  satisfies  a  con¬ 
tracting  recursive  identity  similar  to  (6),  and  is  evaluated  iteratively  in 
the  same  manner  as  p.  Since  this  function  is  given  explicitly  as  part  of 
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the  algorithm  in  the  appendix,  we  avoid  repeating  it  here.  , 

This  spacefilling  curve  satisfies  (PI)  and  (P2)  as  we  now  show.  Note 
that  (6)  breaks  the  range  of  9  into  quarters  to  determine  the  first  bit  of 
the. binary  expansion  of  x  and  y.  Each  successive  pair  of  bits  in  0  corresponds 
to  a  new  bit  for  x  and  a  new  bit  for  y-  Thus,  it  is  cldar  that  satisfies  (PI). 
For  any  9,9',  we  may  construct  a  square  6i  side  at  most  4/j  9—0 * |  containing 

and  so  f(4)*  4/2A,  0  <_  L  h,  in  (P2),  where  {|»|j  in  (3)  is  under¬ 

stood  to  be  the  Euclidean  distance.  We  may  improve  thi3  bound  to 

f  (A)  *  2/T  0  <  A  <  h  (7) 

The  justification  for  (7)  is  too  long  to  be  given  here. 

This  definition  of  $  is  readily  extended  to  d-dimensional  space.  Each 
d  bits  of  0  then  determine  a  single  bit  for  each  coordinate,  and  f  (A)  ■  4  •  S5  * 

If  the  rectilinear  metric  is  to  be  considered,  f  (a)  »  4  »  d  •  for  the 

sup  norm  metric,  f(d)  *  4 


5.  Performance  Analysis 

Computation  Effort.  Since  our  algorithm  projects  the  given  points  onto  C  in  O(kN) 
operations  (in  view  of  (PI)),  and  sorts  them  in  O(NlogN)  operations,  it  requires 
0,(NlogN)  operations  to  obtain  the  heuristic  tour. 

Coding.'  A  short  BASIC  code,  given  in  the  appendix,  demonstrates  the 
ease  with  which  this  algorithm  may  be  implemented. 

Worst-Case  Analysis.  Given  N  sorted  points  the  tour  length  is 

bounded  above  by  ^  f (9^+j-8^) .  Since  this  expression  is  concave 

in  0^,...,®^,  it  achieves  a  maximum  of  Nf(l/N).  In  view  of  (7),  the  heuristic 
tour  cannot  exceed  l/S  in  length.  Projected  onto  a  square  of  area  A,  this  implies 

(8) 


heuristic  tour  <  2/HA 
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An  interesting  corollary  1st  optimal,  tour  <  2/NA. 

Probabilistic  Analysis.  If  the  points  are  uniformly  distributed  in  S,  then 
they  will  be  uniformly  distributed  on  C  as  well,  and  so  are  approximated  by 
a  Poisson  process.  Consequently,  the  expected  tour  length  is  bounded  above 
by  N  J  Ne  f(x)dx  *  An.  Because  of  the  recursive  nature  of  the  algorithm, 
it  is  easy  to  show  that,  in  a  square  of  area  A, 


heuristic  tour 


*>  a  constant,  as  N  4  « 


(9) 


in  much  the  same  sense  that  (optimal  tour  /  KA )  ,765;  see  Beardwood,  Halton 

and  Hammersley[l] .  We  have  estimated  the  parameter  in  (9)  be  .956,  and  so 
the  heuristic  tour  will  be  roughly  25%  over  optimum  when  N  is  large.  Interestingly 
the  ratio  of  heuristic  tour  to  optimal  tour  does  not  depend  on  the  points’ 
distribution  -  so  long  aa  it  has  bounded  density .  We  have  also  shown  that  the 
longest  distance  between  successive  points  in  the  heuristic  tour  is  bounded 
above  by  2/(A/N)  logN ,  almost  surely,  as  N  -+  «;  it  may  generally  be  found 
between  1.1  /(A/N)logN  and  1.3/(A/N)  logN. 

Optional  Enhancement »  An  O(kNlogN)  additional  step  can  reduce  the  expected 
ratio  of  heuristic  tour  to  optimal  tour  to  an  estimated  1.15.  Briefly,  the 
enhancement  projects  each  point  onto  the  boundary  of  its  minimal  containing 
region  and  attempts  to  interchange  its  position  on  C  with  that  of  the  inverse 
image  of  the  point  on  the  opposite  side  of  this  boundary. 


6.  Conclusions 

Our  algorithm  should  prove  useful  in  large  applications  because  of  its 
speed, only  0(N  log  N),  an  order  of  magnitude  faster  than  any  other  TSP 
heuristic  commonly  considered.  It  achieves  this  by  the  surprising  tactic 
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of  ignoring  interpoint  distances  (there  areOCN*')) ,  Its  simplicity  should 
make  it  attractive  for  small  problems  as  well.  The  heuristic,  tour  it  produces 
has  nice  properties:  points  may  be  easily  inserted  or  deleted  (in  O(logN) 
operations)  without  re-solving  the  entire  problem,  and  the  longest  distance 
between  successive  points  on  the  tour  tends  to  be  smal-  —  (This  Is  not  true* 
in  general,  for  tours  that  are  nearest  neighbor,  optimal ,  etc.)  Moreover 
the  tour  tends  to  be  good  with  respect  to  a  variety  of  metrics  and  for 
points  drawn  from  general  (possibly  unknown)  distributions.  The  algorithm 
requires  no  real  multiplications  or  square  roots,  and  so  should  execute 
quickly  on  microprocessor-based  systems.  The  authors  have  even  become 
relatively  proficient  at  solving  small  (up  to  100  points)  problems  by 
hand,  graphically.  Alternatively,  the  points  and  their  images  9  may  be 
marked  on  index  cards  and  placed  in  a  shoebox.  To  produce  a  heuristic 
tour  of  any  subset  of  the  points,  locate  their  cards  and  perform  a  manual 


so  rt. 
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APPENDIX 

The  following  BASIC  program  will  compute  a  heuristic  tour  of  N  points 
in  a  square  of  aide  S.  The  following  remarks  will  help  to  Interpret  it. 
Lines  50-100  compute  the  inverse  image  TH(I)  cf  the  point  (X(I)/S,Y(I)/S) 
in  the  unit  square,  taking  into  account  only  the  first  K*10  bits  of  the 
coordinates1  binary  expansion.  In  the  J-*th  iteration  of  line  70,  the 
lQ(J)-th  quadrant  of  the  present  subsquare  is  selected  as  the  next  sub¬ 
square  to  be  examined.  Lines  110-220  perform  a  sort  to  identify  the  I-th 
smallest  value  in  TH  as  the  ID(I)-th  element  of  the  array.  Lines  150-180 
will  be  executed  N  flog  n]  times.  Variables  whose  names  start  with  I-N 
will  contain  integer  values.  This  program  may  require  simple  modifications 
to  run  under  some  versions  of  BASIC  (e.g. ,  eliminate  variable  dimensions 


in  30). 
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10  INPUT  "NUMBER  OF  POINTS"  ;N  :  K=j.l>  :  KP=2t(K-l) 

20  INPUT  "SIDE  OF  SQUARE" ;S  :  PRINT  "NAME.  X,  Y:" 

30  DIM  &$(N),X(N>,Y(N),TH<N),IQ(IO,1B(N',IO(N) 

40  FOR  1=1  TO  N  :  INPUT  A$(I) ,X(I) ,Y (I)  :  NEXT  I 

50  FOR  1=1  TO  N  :  KR=2*KP-.l 

60  KX= INT (X ( I) *KR/S )  :  KY=INT(Y(I)*KR/S) 

70  FOR  J*1  TO  K  :  JX=INT(KX/KP)  :  JY=INT(KY/KP) 

:  KX=2*(KX-KP*JX)  :  KY=2*(KY-KP*JY) 

:  IQ ( J) = JY+3*JX-2 *JX* JY  :  NEXT  J 
80  T=IQ(K) / 4 

90  FOR  J=K-1  TO  1  STEP  -1  :  T*T+(6-IQ(J))/4 

:  T=T-INT(T)  :  T»(3.5+T+IQ(J))/4  :  NEXT  J 
100  TH(I)=T-INT(T)  :  ID(I)=I  :  NEXT  I 
110  M2=2  :  TH(0)=-1 

120  M=M2  :  M2=M2*2  :  IF  M2<N  THEN  GOTO  120 
130  FOR  I=N  TO  N-M+l  STEP  -1 
140  J1=I  :  J2*I-M  :  J3*0  :  J4=I 
150  J3=J3+1 

160  IF  J1<0  THEN  J1=0  ELSE  IF  J2<0  THEN  J2=0 
170  IF  TH(ID(J1) ) >TH(ID( J2) )  THEN  IC(J3)=ID(J1) 

:  J1-J1-M2  ELSE  IC(J3)=ID<J2)  :  J2=J2-M2 
180  IF  J1>0  OR  J2>0  THEN  GOTO  150 

190  IF  J3<»1  THEN  GOTO  220 

200  ’  FOR  J=1  TO  J3  :  ID(J4T=IC(J)  :  J4*=J4-M  :  NEXT  J 

210  NEXT  I 

220  M2-M  :  M«M/2  :  IF  M>=1  THEN  GOTO  130 

230  PRINT  :  PRINT  "RANK" , "NAME” , "X" , "Y" , "THETA" 

240  FOR  I«1  TO  N  :  J«ID(I) 

'  :  PRINT  I,A$(J),X(J),Y(J),TH(J)  :  NEXT  I 


250  END 
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Figure  1:  Recursive  Structure  of  the  Spacefilling 
Curve  Given  by  (6) 


Figure  2:  Successive  Approximations  to  i (i 


